% quadrature 
[Cheb_nodes,~] = qnwcheb(nE,-1,1);
a           = -5;
b           = 5;
Eshocks     = (b-a)*(Cheb_nodes'+1)/2 + a;
temp        = pi*(b-a)/(2*nE);
weights     = temp.*sqrt(1-Cheb_nodes'.^2);
probE       = normpdf(Eshocks).*weights;
probE           = probE';
Eshocks         = Eshocks';


% ATM asset put
KS_value = 1;
lamAx = lamA - (1-tau_c)*(1-tau_d); 
g = exp(mu_XQ + sig_X.*Eshocks');
ATM_Asset_Put = NaN(nS,1);
ATM_Asset_IV = NaN(nS,1);
for is=1:nS % today's state
    ATM_Asset_Put(is) = (Q(is,:) * max(0,lamAx(is)*KS_value - lamAx*g(is,:)) * probE)/Rf(is,1);
    PV_div = (g(is,:)* probE)/Rf(is,1);
    ATM_Asset_IV(is) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_value*lamAx(is),log(Rf(is,1))*12,1/12,vol)-ATM_Asset_Put(is),0.2);
end


% IV skew in ATM asset IV units
nKS = 31;
MON_values = linspace(-2,1,nKS)';
KS_values = exp(MON_values * (ATM_Asset_IV'/sqrt(12)));
Asset_Put = NaN(nS,nKS);
Asset_IV = NaN(nS,nKS);
for is=1:nS % today's state
    PV_div = (g(is,:)* probE)/Rf(is,1);
    for iks=1:nKS
        Asset_Put(is,iks) = (Q(is,:) * max(0,lamAx(is)*KS_values(iks,is) - lamAx*g(is,:)) * probE)/Rf(is,1);
        Asset_IV(is,iks) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_values(iks,is)*lamAx(is),log(Rf(is,1))*12,1/12,vol)-Asset_Put(is,iks),0.2);
    end
end


% Asset-to-earnings ratio
% eps_it : Idiosyncratic shock, corr with eps_ct is rho (under P and Q)
X           = (exp(mu_XQ + sig_XS.*Eshocks')*probE) ./ Rf(:,1); % [nS,1]
lamA_noIdio = (1-tau_c)*(1-tau_d)*((eye(nS)-X.*Q)\ones(nS,1)); % [nS,1]


% ATM asset put
KS_value = 1;
lamAx = lamA_noIdio - (1-tau_c)*(1-tau_d); 
g = exp(mu_XQ + sig_XS.*Eshocks');
ATM_Asset_Put_noIdio = NaN(nS,1);
ATM_Asset_IV_noIdio = NaN(nS,1);
for is=1:nS % today's state
    ATM_Asset_Put_noIdio(is) = (Q(is,:) * max(0,lamAx(is)*KS_value - lamAx*g(is,:)) * probE)/Rf(is,1);
    PV_div = (g(is,:)* probE)/Rf(is,1);
    ATM_Asset_IV_noIdio(is) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_value*lamAx(is),log(Rf(is,1))*12,1/12,vol)-ATM_Asset_Put_noIdio(is),0.2);
end


% IV skew in ATM asset IV units
nKS = 31;
MON_values = linspace(-2,1,nKS)';
KS_values = exp(MON_values * (ATM_Asset_IV_noIdio'/sqrt(12)));
Asset_Put_noIdio = NaN(nS,nKS);
Asset_IV_noIdio = NaN(nS,nKS);
for is=1:nS % today's state
    PV_div = (g(is,:)* probE)/Rf(is,1);
    for iks=1:nKS
        Asset_Put_noIdio(is,iks) = (Q(is,:) * max(0,lamAx(is)*KS_values(iks,is) - lamAx*g(is,:)) * probE)/Rf(is,1);
        Asset_IV_noIdio(is,iks) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_values(iks,is)*lamAx(is),log(Rf(is,1))*12,1/12,vol)-Asset_Put_noIdio(is,iks),0.2);
    end
end

